function [SVAN,SIGMA]=swstate(S,T,P0,ftn);

% SWSTATE State equation for seawater

%

%        [SVAN,SIGMA]=SWSTATE(S,T,P) returns the specific volume

%        anomaly SVAN (m^3/kg*1e-8) and the density anomaly SIGMA (kg/m^3)

%        given the salinity S (ppt), temperature T (deg C) and pressure

%        P (dbars).

%

%        [dVdT,dRdT]=SWSTATE(S,T,P,'dT') returns derivatives w.r.t.

%        temperature of the volume and density.

%

%        [dVdS,dRdS]=SWSTATE(S,T,P,'dS') returns derivatives w.r.t.

%        salinity.

%

%        [dVdP,dRdP]=SWSTATE(S,T,P,'dP') returns derivatives w.r.t.

%        pressure.

%

%        All elements can be scalars, vectors, or matrices but should be

%        the same size.



%Notes: RP (WHOI 2/dec/91)

%

%  This stuff is directly copied from the UNESCO algorithms, with some

%  minor changes to make it Matlab compatible (like adding ";" and changing

%  "*" to ".*" when necessary.

%

%      RP  (WHOI 3/dec/91)

%

%  Added first derivative calculations.



derivT=0;

derivS=0;

derivP=0;



if (nargin==4),

   if     (ftn=='dT'), derivT=1;

   elseif (ftn=='dS'), derivS=1;

   elseif (ftn=='dP'), derivP=1;

   else error('swstate: Unrecognized option!');

   end;

end;



% ******************************************************

% SPECIFIC VOLUME ANOMALY (STERIC ANOMALY) BASED ON 1980 EQUATION

% OF STATE FOR SEAWATER AND 1978 PRACTICAL SALINITY SCALE.

% REFERENCES

% MILLERO, ET AL (1980) DEEP-SEA RES.,27A,255-264

% MILLERO AND POISSON 1981,DEEP-SEA RES.,28A PP 625-629.

% BOTH ABOVE REFERENCES ARE ALSO FOUND IN UNESCO REPORT 38 (1981)

% UNITS:      

%       PRESSURE        P0       DECIBARS

%       TEMPERATURE     T        DEG CELSIUS (IPTS-68)

%       SALINITY        S        (IPSS-78)

%       SPEC. VOL. ANA. SVAN     M**3/KG *1.0E-8

%       DENSITY ANA.    SIGMA    KG/M**3

% ******************************************************************

% CHECK VALUE: SVAN=981.3021 E-8 M**3/KG.  FOR S = 40 (IPSS-78) ,

% T = 40 DEG C, P0= 10000 DECIBARS.

% CHECK VALUE: SIGMA = 59.82037  KG/M**3 FOR S = 40 (IPSS-78) ,

% T = 40 DEG C, P0= 10000 DECIBARS.

%HECK VALUE: FOR S = 40 (IPSS-78) , T = 40 DEG C, P0= 10000 DECIBARS.

%        DR/DP                  DR/DT                 DR/DS

%       DRV(1,7)              DRV(2,3)             DRV(1,8)

%

% FINITE DIFFERENCE WITH 3RD ORDER CORRECTION DONE IN DOUBLE PRECSION

%

%       3.46969238E-3       -.43311722           .705110777

%

% EXPLICIT DIFFERENTIATION SINGLE PRECISION FORMULATION EOS80 

% 

%       3.4696929E-3        -.4331173            .7051107

%

% (RP...I think this ---------^^^^^^ should be -.4431173!);





% *******************************************************

% DATA

   R3500=1028.1063;

   R4=4.8314E-4;

   DR350=28.106331;



% CONVERT PRESSURE TO BARS AND TAKE SQUARE ROOT SALINITY.

      P=P0/10.;

      SAL=S;

      SR = sqrt(abs(S));

% *********************************************************

% PURE WATER DENSITY AT ATMOSPHERIC PRESSURE

%   BIGG P.H.,(1967) BR. J. APPLIED PHYSICS 8 PP 521-537.

%

      R1 = ((((6.536332E-9*T-1.120083E-6).*T+1.001685E-4).*T  -9.095290E-3).*T+6.793952E-2).*T-28.263737;

% SEAWATER DENSITY ATM PRESS. 

%  COEFFICIENTS INVOLVING SALINITY

      R2 = (((5.3875E-9*T-8.2467E-7).*T+7.6438E-5).*T-4.0899E-3).*T+8.24493E-1; 

      R3 = (-1.6546E-6*T+1.0227E-4).*T-5.72466E-3;

%  INTERNATIONAL ONE-ATMOSPHERE EQUATION OF STATE OF SEAWATER

      SIG = (R4*S + R3.*SR + R2).*S + R1;

% SPECIFIC VOLUME AT ATMOSPHERIC PRESSURE

      V350P = 1.0/R3500;

      SVA = -SIG*V350P./(R3500+SIG);

      SIGMA=SIG+DR350;

      V0 = 1.0./(1000.0 + SIGMA);

%  SCALE SPECIFIC VOL. ANAMOLY TO NORMALLY REPORTED UNITS

      SVAN=SVA*1.0E+8;



if (derivS),               % These are derivatives for (S,T,0).

      R4S=9.6628E-4;

      RHO1 = 1000.0 + SIGMA;



      RHOS=R4S*SAL+1.5.*R3.*SR+R2;  

      V0S=-RHOS./(RHO1.*RHO1);  

elseif (derivT),

      R1 =(((3.268166E-8*T-4.480332E-6).*T+3.005055E-4).*T -1.819058E-2).*T+6.793952E-2;

      R2 = ((2.155E-8*T-2.47401E-6).*T+1.52876E-4).*T-4.0899E-3;

      R3 = -3.3092E-6*T+1.0227E-4;

      RHO1 = 1000.0 + SIGMA;



      RHOT = (R3.*SR + R2).*SAL + R1;    

      V0T = -RHOT./(RHO1.*RHO1);

end;

      

% ******************************************************************

% ******  NEW HIGH PRESSURE EQUATION OF STATE FOR SEAWATER ********

% ******************************************************************

%        MILLERO, ET AL , 1980 DSR 27A, PP 255-264

%               CONSTANT NOTATION FOLLOWS ARTICLE

%********************************************************

% COMPUTE COMPRESSION TERMS

      E = (9.1697E-10*T+2.0816E-8).*T-9.9348E-7;

      BW = (5.2787E-8*T-6.12293E-6).*T+3.47718E-5;

      B = BW + E.*S;    % Bulk Modulus (almost)

%  CORRECT B FOR ANAMOLY BIAS CHANGE

      Bout = B + 5.03217E-5;



if (derivS),

      DBDS=E;

elseif (derivT),

      BW = 1.05574E-7*T-6.12293E-6;

      E = 1.83394E-9*T +2.0816E-8;

      BT = BW + E.*SAL;

end;

%             

      D = 1.91075E-4;

      C = (-1.6078E-6*T-1.0981E-5).*T+2.2838E-3;

      AW = ((-5.77905E-7*T+1.16092E-4).*T+1.43713E-3).*T-0.1194975;

      A = (D*SR + C).*S + AW;    

%  CORRECT A FOR ANAMOLY BIAS CHANGE

      Aout = A + 3.3594055;



if (derivS),

      DADS=2.866125E-4*SR+C;

elseif (derivT),

      C = -3.2156E-6*T -1.0981E-5;

      AW = (-1.733715E-6*T+2.32184E-4).*T+1.43713E-3;

      AT = C.*SAL + AW;

end;

           

      B1 = (-5.3009E-4*T+1.6483E-2).*T+7.944E-2;

      A1 = ((-6.1670E-5*T+1.09987E-2).*T-0.603459).*T+54.6746;

      KW = (((-5.155288E-5*T+1.360477E-2).*T-2.327105).*T+148.4206).*T-1930.06;

      K0 = (B1.*SR + A1).*S + KW;



if (derivS),

      K0S=1.5*B1.*SR+A1;

      KS=(DBDS.*P+DADS).*P+K0S;

elseif (derivT),

      B1 = -1.06018E-3*T+1.6483E-2;

      % APRIL 9 1984 CORRECT A1 BIAS FROM -.603457 !!!

      A1 = (-1.8501E-4*T+2.19974E-2).*T-0.603459;

      KW = ((-2.0621152E-4*T+4.081431E-2).*T-4.65421).*T+148.4206;

      K0T = (B1.*SR+A1).*SAL + KW;

      KT = (BT.*P + AT).*P + K0T;

end;





% EVALUATE PRESSURE POLYNOMIAL 

% ***********************************************

%   K EQUALS THE SECANT BULK MODULUS OF SEAWATER

%   DK=K(S,T,P)-K(35,0,P)

%  K35=K(35,0,P)

% ***********************************************

      DK = (B.*P + A).*P + K0;

      K35  = (5.03217E-5*P+3.359406).*P+21582.27;

      GAM=P./K35;

      PK = 1.0 - GAM;

      SVA = SVA.*PK + (V350P+SVA).*P.*DK./(K35.*(K35+DK));

%  SCALE SPECIFIC VOL. ANAMOLY TO NORMALLY REPORTED UNITS

      SVAN=SVA*1.0E+8;      % Volume anomaly

      V350P = V350P.*PK;

%  ****************************************************

% COMPUTE DENSITY ANAMOLY WITH RESPECT TO 1000.0 KG/M**3

%  1) DR350: DENSITY ANAMOLY AT 35 (IPSS-78), 0 DEG. C AND 0 DECIBARS

%  2) DR35P: DENSITY ANAMOLY 35 (IPSS-78), 0 DEG. C ,  PRES. VARIATION

%  3) DVAN : DENSITY ANAMOLY VARIATIONS INVOLVING SPECFIC VOL. ANAMOLY

% ********************************************************************

% CHECK VALUE: SIGMA = 59.82037  KG/M**3 FOR S = 40 (IPSS-78),

% T = 40 DEG C, P0= 10000 DECIBARS.

% *******************************************************

      DR35P=GAM./V350P;

      DVAN=SVA./(V350P.*(V350P+SVA));

      SIGMA=DR350+DR35P-DVAN;  % Density anomaly



      K=K35+DK;

      VP=1.0-P./K;

      V = (1.) ./(SIGMA+1000.0);



if (derivS),

      VS=V0S.*VP+V0.*P.*KS./(K.*K);



      SVAN=VS;              % dVdS

      SIGMA=-VS./(V.*V);    % dRdS

elseif (derivT),

      VT = V0T.*VP + V0.*P.*KT./(K.*K);



      SVAN=VT;              % dVdT

      SIGMA=-VT./(V.*V);    % dRdT

elseif (derivP),

      DKDP = 2.0*Bout.*P + Aout;   

% CORRECT DVDP TO PER DECIBAR BY MULTIPLE *.1

      DVDP = -.1*V0.*(1.0 - P.*DKDP./K)./K;



      SVAN=DVDP;            % dVdP

      SIGMA=-DVDP./(V.*V);  % dRdP

end;









